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Abstract. The evaluation of a matrix exponential function is a classic problem of 
computational linear algebra. Many different methods have been employed for its 
numerical evaluation [Moler C and van Loan C 1978 SI AM Review 20 4], none of which 
produce a definitive algorithm which is broadly applicable and sufficiently accurate, as 
well as being reasonably fast. Herein, we employ a method which evaulates a matrix 
exponential as the solution to a first-order initial value problem in a fictitious time 
variable. The new aspect of the present implementation of this method is to use finite 
elements in the fictitious time variable. [Weatherford C A, Red E, and Wynn A 2002 
Journal of Molecular Structure 592 47] Then using an expansion in a properly chosen 
time basis, we are able to make accurate calculations of the exponential of any given 
matrix as the solution to a set of simultaneous equations. 



Introduction 

The evaluation of the exponential of a square matrix e'^ is a classic problem of 
computational linear algebra. [1] A large number of methods have been proposed and 
used for its evaluation. None of these methods have produced a generally applicable 
method which is sufficiently accurate as well as being reasonably fast. Thus this problem 
might be considered unsolved. It is clearly an important and pervasive problem which 
arises in a wide variety of contexts. [2,3] A method which is capable of producing a 
solution to essentially arbitrary precision would thus be of great importance. The 
present work uses a variant of a method described in Ref. 1, namely the introduction of 
an artificial time-parameter which produces an initial- value problem. Instead of calling 
a 'canned' solver, the present work uses a method introduced by several of the present 
authors [4] to solve quantum mechanical initial-value problems. In this method, a finite 
element technique is used to propagate from an initial condition at t = 0, which is the 
unit matrix, to the desired result at t = 1. The time axis is broken up into an arbitrary 
number of time elements and the solution is propagated from element to element, using 
a special basis in time introduced here for the first time. 
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The next section presents the analysis of the problem and describes the solution 
algorithm. Then the algorithm is applied to evaluate the exponential of a number of 
test matrices. Finally, the conclusions are presented. 



Analysis and Solution Algorithm 

The problem at hand is the evaluation of the exponential of a square (generally complex) 
matrix e^. The present method introduces an artificial time parameter so as to 
transform the evaluation into the solution of an initial-value problem. For a given 
n X n square matrix A, consider the following parametrized function definition: 

m ^ e^*, (1) 

where ^ is also an n x n square matrix and the desired solution is ^'(1) = which 
evolves from the initial value given by ^'(0) = 1 (1 is the diagonal unit matrix). This is 
a solution of the following linear ordinary differential equation, written for each matrix 
element, 

n 

^ij{t) = YlAk'^kj{t)^ (2) 

k=l 

where the over-dot stands for the time-derivative. The time axis extends over the interval 
[0, 1]. Now break the time axis up into elements that extend between nodes ti and ii+i, 
and define a local time r that spans [—1,1]. The local time transformation is defined 
by the relation, 

T^qt-p, (3) 

where, q = 2/(ti_|_i —ti) and p = {ti+i+ti) / (ti+i — ti). Thus, for an arbitrary time element 
e, Eq. (2) can be written in terms of local time r as 

QMf{r) = tA^Mhr). (4) 

k=l 

At this point, we will use the following ansatz for to enforce continuity between 
two consecutive finite elements 

*g^(r) = fifir) + *g-^^(+l), iJH-l) = (5) 



itj 

m—l 



and expand fifir) as 

f^fir) = E B^^^^^s.ir) (6) 
in a basis we define by 



f_T,{T)d.T (7) 
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where T^j^ij) are Chebyshev Polynomials of the first kind. [5] Note that these basis 
functions enforce the initial condition on the /'s given in Eq. (5) since 1) — 0. 
The result for the decomposition of / in m basis functions is 

m— 1 
m—l 

= E B'^^'Kir)- (9) 

Now, insert (8) and (9) into Eq. (4), and multiply from the left by w{r)s^/{T) and 
integrate from —1 to +1 (note that w{t) = (1 — r^)"-'^/^ is the weighting function for 
Chebyshev polynomials). Rearranging terms we get, 



or, rearrangmg 



(10) 



?E[/' V(^VWr.Mdr]Bf) =^A;^)[/' s,,(r)a;(r)s,(r)dr]BjKe) 

+ E 4! [ /' V irMr)To{r) drj^g"^) (+1) 
where, 7o(r) = 1. Defining, the integrals in the above equation as 

C,,,^ f\,,iTMr)T,{r)dr (11) 
= J'^ s^,{t)u;{t)s^{t) dr (12) 

V = J_^Si^>{T)uj{T)To{T)dT (13) 

and substituting Eqs. (11-13) into Eq. (10) gives, 

<1 E C,',B^^ = E 4^D,,,B^^^ + V E t (+1) (14) 

IJ- kn k 



fj,k k 

where 6ik is the usual Kronecker delta function. Then rewrite Eq. (15) as 
where 



(17) 
(18) 
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Equation (16) is a set of simultaneous equations of size {n x m), which can be 
written in matrix form as, 

Q{e)-Qi(^) ^ pj(e,e-l) j ^ ^ 3, 71. (19) 

Here, f2*^^) is a (complex) matrix and for each j, T^i^'^-'^^ and B-"^^) are vectors. Eq. (19) 
applies for each time element e. The solution is propagated from element to element, 
from t — Otot — 1. The above equation can be solved numerically in many ways, but we 
have chosen the method of LU decomposition. [6] The present method is ideally suited 
to high-performance computers where the solver of choice would probably be iterative. 
In the present case, we apply LU decomposition to fl'^'^^ and back substituting all of 
the r'^('^''^~i)'s, we will have all the elements for the matrix (which can also be viewed as 
three dimensional) This LU decomposition only needs to be done once since fl^'^^ 

is independent of time. Thus, the propagation just involves a matrix vector multiply. 
Then, we employ Eq. (8) to solve for ^(^^(r = 1) for the element e, which, in turn, will 
be used as ^(^+i)(r = —1) for the next element e + 1. Starting off with a unit matrix 
for "if^^^t = 0), we continue this process till we calculate *(t = 1) at the last node, 
which is the exponential of the given matrix A. 

Results 

The calculations presented below were done on a Macintosh Intel laptop using Gnu C++ 
which has machine accuracy limit of 2.22045 x 10~^^. As an illustration, let's borrow a 
'pathological' matrix from [1], which we have modified slightly to make it even worse. 
Consider a matrix Ml given by. 



Ml = 



-73 
-96 



36 
47 





" -1 " 




"13" 




-25 




2 4 



1 3 

2 4 

The exponent of Ml can be easily calculated as, 

1 3 

2 4^^ 

-1 



(20) 



Ml 





" ■ 




" -2 3/2 




_ 




1 -1/2 



— e 

2e-25 



25^1 



(21) 



-2e-^ + 3e-25 (3/2) (e" 
-4e-i + 4e-25 Se'^ - 
The above matrix, exact to 16 decimal places, is given by 

-0.7357588823012208 0.5518191617363316 
-1.4715177646302175 1.1036383234865511 

The result of our program is displayed below and we run it by using just 8 time steps 
and 8 basis functions. The result is accurate to 13 decimal places already. 

-0.7357588823012(181) 0.5518191617363(358) 
-1.4715177646302(120) 1.1036383234865(592) 



gMl 



(22) 
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As an example of a non - diagonalizable matrix, consider the following matrix M2, 
with complex eigenvalues 



M2 



It can be shown that, 



-1 

1 



cos(l) 
sin(l) 



—sin{l] 
cos(l) 



(23) 



(24) 



We are able to achieve 14 decimal digit accuracy with 8 time steps and 8 basis functions. 

(25) 



gM2 



0.54030230586814 -0.84147098480790 
0.84147098480790 0.54030230586814 



Table 1 shows the minimum number of basis functions, for a given number of time steps, 
which were required to achieve a precision of ±1 x 10~^^ on matrices whose exponential 
is known exactly. The matrices chosen are: the simplest possible matrix - a 2 x 2 real 
unit matrix, for which the result is the constant e on the diagonals, and the matrices 
Ml and M2. 

Let's check our program on matrices, which wc picked randomly and for which we 
had no apriori knowledge as to the result of their exponentiation. Wc fixed 8 time steps 
and/or 8 basis functions, and varied the other corresponding parameter from 5 to 40 
and checked how the results of the program varied in accuracy. For the sake of saving 
space, we only displayed the result of the last element-the other elements of the matrix 
exponential behaved similarly. The matrices chosen are a 5 x 5 real matrix M3, 



M3 = 



and a 3 X 3 complex matrix M4 



-0.1 


-0.2 


-0.3 


-0.4 


-0.5 


-0.6 


-0.7 


-0.8 


-0.9 


-1 


0.1 


0.2 


0.3 


0.4 


0.5 


0.6 


0.7 


0.8 


0.9 


1 


1 


2 


3 


4 






(26) 



M4 



1 + i 1 — i i 
1 2i 
l + 2i + i 



(27) 



From Table 2, one can see that the numbers up to 12 decimal digits have saturated 
after 5 time steps and/or basis functions. Similarly, Table 3 shows 13 digits of accuracy 
as we switch the two parameters from 5 to 40, except for the case of 5 basis functions, 
which only shows 8 accurate significant digits. This shows that for complex matrices, 
there is inherently more work for the program to handle because of the imaginary part 
of the matrix elements and there is apparently more sensitivity to the number of basis 
functions used than to the number of time steps. 
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Table 1. Minimum number of basis functions and time steps required for a precision 
of ±1 X 10"" for a 2 X 2 unit matrix, Ml and M2. 



2x2 unit matrix Ml M2 

Time steps Basis functions Time steps Basis functions Time steps Basis functions 



1 


11 


05 




1 


11 


2 


9 


08 


7 


2 


9 


4 


8 


16 


6 


4 


8 


8 


7 


50 


5 


8 


7 


16 


6 


256 


4 


15 


6 


58 


5 






40 


5 



Table 2. Results of matrix e^^ for typical runs of 8 time steps and 8 basis functions. 



Result 


Time steps 


Basis functions 


3.210309305973118 


5 


8 


3.210309305973288 


40 


8 


3.210309315373377 


8 


5 


3.210309305973281 


8 


40 



Table 3. Results of matrix 633 for typical runs of 8 time steps and 8 basis functions. 



Result 




Time steps 


Basis functions 


-0.511977122298063 - 


i 0.089772811313512 


5 


8 


-0.511977122298081 - 


i 0.089772811313526 


40 


8 


-0.511977121264660- 


i 0.089772810979965 


8 


5 


-0.511977122298082 - 


i 0.089772811313526 


8 


40 



Conclusion 

Wc have presented a robust, easily used, and accurate algorithm for the evaluation of 
the exponential of a matrix. We did this by introducing an artificial time parameter 
and evaluating the matrix exponential as the solution of an initial-value problem in this 
artificial time. We solved the initial-value problem by using finite elements in time with 
a new time basis which we defined here so as to enforce the initial conditions on the so- 
lution at the beginning of each time finite element. This resulted in set of simultaneous 
equations for the expansion coefficients. The actual algorithm employed here was an LU 
decomposition which was very fast and efficient. The relative efficiency of the method 
should be most apparent when implemented on high-performance computers since the 
algorithm is highly parallel. The method was applied to several matrices as a proof of 
the validity of the algorithm. The results of our calculations show that we only need 
about 8 basis functions and 8 time steps for the matrices considered for accuracies as 
great as 13 significant digits. We trust that this method of numerically calculating the 
exponential of a matrix will be recognized to be a nondubious one! 
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